Fourier's Law in a Quantum Spin Chain and the Onset of Quantum Chaos 
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We study heat transport in a nonequilibrium steady state of a quantum interacting spin chain. We 
provide clear numerical evidence of the validity of Fourier law. The regime of normal conductivity 
is shown to set in at the transition to quantum chaos. 
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The derivation of Fourier's law of heat conduction from 
the microscopic dynamics, without any ad hoc statistical 
assumption, is one of the great challenges of nonequilib- 
rium statistical mechanics Even in the context of 
classical dynamical systems, the issue of energy (heat) 
transport, in spite of having a long history (recently re- 
viewed in 0]), is not completely settled. It has become 
clear that transport theory requires that the underly- 
ing deterministic dynamics yield a truly diffusive pro- 
cess. On the other hand, it is known that classical non- 
linear systems of interacting particles, above a critical 
interaction strength, typically exhibit chaotic and diffu- 
sive behaviour which then leads to the onset of Fourier 
law J — —kVT, relating the macroscopic heat flux to 
the temperature gradient VT. Therefore deterministic 
chaos appears to be an essential ingredient required by 
transport theory. In this perspective, of particular in- 
terest is the problem, almost completely unexplored, of 
the derivation of Fourier law from quantum dynamics. 
So far, investigations have been mainly focused on lin- 
ear response theory)^ la L3, The possibility instead 
to derive the Fourier law from quantum dynamics, by 
establishing the dependence of J on VT in a nonequi- 
librium steady state, calls directly in question the issue 
of quantum chaos. In this connection, a main feature 
of quantum motion is the lack of exponential dynamical 
instability 9] , a property which is at the heart of classical 
dynamical chaos. This fact may render very questionable 
the possibility to derive the Fourier law of heat conduc- 
tion in quantum mechanics. However, quite interestingly, 
it has been shown that strong, exponential unstable, clas- 
sical chaos is not necessary |2( (actually, strictly speaking, 
is not even sufficient for normal transport. Thus it 
is interesting to inquire if, and under what conditions, 
Fourier law emerges from the laws of quantum mechan- 
ics. This is the purpose of the present paper. 

To investigate this problem one has to deal with a finite 
open system connected to heat baths. Here we consider 
an interacting quantum spin- 1/2 chain which exhibits the 
transition from integrability to quantum chaos as a pa- 
rameter, e.g. the magnetic field, is varied. The stan- 
dard treatment of this problem is based on the master 



equation, thus limiting investigations to relatively small 
system sizes. By using this method, in an interesting 
paper |10|. the decay of current correlation function in a 
model of non-integrable chain of quantum spins is com- 
puted. Here we take a different approach namely we fol- 
low the evolution of the system described by a pure state, 
which is stochastically coupled to an idealized model of 
heat baths. Stochastic coupling is realized in terms of 
a local measurement at the boundary of the system and 
stochastic but unitary exchange of energy between the 
system and the bath. By this method we have been able 
to perform very effective numerical simulations which al- 
low to observe a clear energy /temperature profile and to 
measure the heat current J. In the nonintegrable regime 
where the spectral statistics is described by random ma- 
trix theory (RMT) - the regime of "quantum chaos" - we 
found very accurate Fourier law scaling J /AT oc 1/L, 
where L is the size of the chain. In the integrable and 
near-integrable regimes instead, we found that the heat 
transport is ballistic J oc L°. 

We consider an Ising chain of L spins 1/2 with cou- 
pling constant Q subject to a uniform magnetic field 
h = (h x , 0, h z ), with open boundaries. The Hamiltonian 
reads 



L-2 

H = -Q a~o 



n=0 



L-l 



(1) 



where the operators a n — (c^, cr^, <r*) are the Pauli ma- 
trices for the n-th spin, n = 0, 1, . . . L — 1. We set the 
coupling constant Q = 2. In this system, the only triv- 
ial symmetry is a reflection symmetry, a n — > c?x,_i_ n . 
Moreover the direction of the magnetic field affects the 
qualitative behavior of the system: If h z =0, the Hamil- 
tonian Q corresponds to the Ising chain in a transversal 
magnetic field. In this case the system is integrable as 
(fT|l can be mapped into a model of free fermions through 
standard Wigner- Jordan transformations. When h z is 
increased from zero, the system is no longer integrable. 
When h z is of the same order of h x quantum chaos sets 
in leading to a very complex structure of quantum states 
as well as to fluctuations in the spectrum that are sta- 
tistically described by RMT 11] . The system becomes 
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FIG. 1: Nearest neighbor level spacing distribution P(s) for 
the chaotic (left), intermediate (center) and integrable (right) 
chains. P(s) was obtained by diagonalizing the Hamiltonian 
([Ql for a chain of length L = 12, averaging over even and odd 
parity subspectra. The curves correspond to Pwd (dashed 
line) and to Pp (dotted line). 



again (nearly) integrable when h z 3> h x . Therefore, by 
choosing the direction of the external field we can explore 
different regimes of quantum dynamics. 

The onset of quantum chaos is commonly studied in 
terms of the nearest neighbor level spacing distribution 
P(s) that gives the probability density to find two ad- 
jacent levels at a distance s. For an integrable system 
the distribution P(s) has typically a Poisson distribution 
Pp(s) = cxp (— s). In contrast, in the quantum chaos 
regime (for Hamiltonians obeying time-reversal invari- 
ance), P(s) is given by the Gaussian Orthogonal En- 
semble of random matrices (GOE). In this case, P{s) 
is well-approximated by the Wigner surmise Pwd(s) — 
(ns/2) exp(— 7rs 2 /4), exhibiting the so-called "level repul- 
sion" . 

In Fig. n w e show the results of our numerical simula- 
tions of system Q for three different values of the mag- 
netic field: (i) chaotic case h = (3.375,0,2) at which the 
distribution P(s) agrees with GOE and thus corresponds 
to the regime of quantum chaos, (ii) integrable case 
h = (3.375,0,0), at which P(s) is close to the Poisson 
distribution, and (Hi) intermediate case h = (7.875,0,2) 
at which the distribution P(s) shows a combination of 
(weak) level repulsion and an exponential tail. We have 
consistently observed signatures for the onset of quantum 
chaos also in other spectral correlations as well as in the 
structure of eigenfunctions. 

Let us now turn to study energy transport in this 
model system. To this end we need to couple both ends 
of the chain of spins to thermal reservoirs at temperature 
T. We have devised a simple way to simulate this cou- 
pling, namely the state of the spin in contact with the 
bath is statistically determined by a Boltzmann distri- 
bution with parameter T . Our model for the reservoirs 
is analogous to the stochastic thermal reservoirs used in 
classical simulations JpJ and we thus call it a quantum 
stochastic reservoir. We use units in which Planck and 
Boltzmann constants are set to unity K = ks = P 



In the representation basis of the wave function at 
time t can be written as 

m))= J2 C-o,.i,...,«-i(*)l»o,*i...ai-i> , (2) 

So,Sl,...,SZ,_l 

where s n = 0, 1 represents the up, down state of the 
rt-th spin, respectively. The wave function at time t 
is obtained from the unitary evolution operator U(i) = 
exp(— iTLt). The interaction with the reservoir is not in- 
cluded in the unitary evolution. Instead, we assume that 
the spin chain and the reservoir interact only at discrete 
times with period r at which the states of the leftmost 
(sq) and the rightmost (s^_i) spins are stochastically re- 
set. Thus, the evolution of the wave function from time 
t to time t + r can be represented as 

Mi + r))=S(A,&)U(r)Mt)) , (3) 

where S(/3i,/3 r ) represents the unitary stochastic action 
of the interaction with the left and right reservoirs at 
temperatures /3j _1 and (3~ x respectively. 

The action of 2(/3i, j3 r ) takes place in several steps: 

(i) The wave function is first rotated by the angle 
a = t&T\T 1 (h x /h z ) to the eigenbasis of the components 
<Ti = h ■ do/h, er r = h ■ <j£_i//i of the edge spins along 
the field, that is \ip) -> e~ iQ(<T ° +CT "-i )/2 |</'). Here, h = \h\ 
stands for the strength of the magnetic field. 

(ii) A local measurement of the observables a\, cr r is per- 
formed. Then the state of the spins at the borders col- 
lapses to a state (sq,s^ j _ 1 ) with probability 

P(4,S*L-l)= E \Cs*, Sl ,..., SL ^,sl_ 1 \ 2 ■ (4) 

Sl,...,S£_2 

In other words this means that we put all coefficients 
C SOtSlt ... tSL _ 1 with (s ,s L _i) ^ (sq- s l-i) to zer0 - 
(Hi) The new state of the edge spins is stochastically cho- 
sen. After this, which simulates the thermal interaction 
with the reservoirs, each of the edge spins is set to down, 
(up) state with propability fj,,(l — /i). The probability 
fi(/3) depends on the canonical temperature of each of 
the thermal reservoirs: 

Mi)= e-frh + ^h 5 ie{l,r}. (5) 

(iv) Finally, the wave function is rotated back to the 
basis, \il>) -> e ia ^ + <-i )/2 \i>). 

This completes the description of the interaction with 
the quantum stochastic bath. This interaction thus (pe- 
riodically) resets the value of the local energy hai jr of the 
spins in contact with the reservoirs. This information 
is then transmitted to the rest of the system during its 
dynamical evolution and relaxation towards equilibirum. 
Therefore, the value of r controls the strength of the 
coupling to the bath. We have found that, in our units, 
r = 1 provides an optimal choice. We have nevertheless 
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FIG. 2: Mean local energy E as a function of the temper- 
ature T for the chaotic chain, for L = 6,8, 10 with crosses, 
triangles and circles, respectively. Each point corresponds to 
the equilibrium simulation in which the bath temperatures 
are set to the same value T. The average was obtained from 
the center of the profile so that, the two sites closest to each 
border were not considered. The dashed curve corresponds 
to the canonical average E ca , n (T) of the local 2-body energy 
operator H n . The two solid squares indicate the values of the 
temperature that were later used in the out of equilibrium 
simulations. In the inset the same data are shown for the 
integrable chain (squares) for L = 8. 



FIG. 3: Energy profile (H n ) for an out of equilibrium situa- 
tion for the chaotic chain. The nominal values of the temper- 
atures of the reservoirs Tj = 5 and T T = 50, are both in the 
high temperature regime (see the solid squares in Fig. [5] The 
different symbols correspond to chains of size L = 8 (crosses) , 
L — 14 (open circles) and L = 20 (solid circles). The dashed 
fine was obtained from a linear fit of the data for L = 20 in 
which the two sites closest to each border were not considered. 
In the insets (I) and (II) the energy profile is shown for the 
integrable and intermediate cases respectively, for L — 15. 



performed simulations for other values of r with qualita- 
tively similar results. 

The unitary evolution U(t)\ip) of the system has been 
computed by an accurate high order split-step factoriza- 
tion of the unitary evolution operator |l3l| . For each run 
the initial wavefunction |^(0)) of the system is chosen at 
random. The system is then evolved for some relaxation 
time T re i after which it is assumed to fluctuate around a 
unique steady state. Measurements are then performed 
as time averages of the expectation value of suitable ob- 
servables. We further average these quantities over dif- 
ferent random realizations of "quantum trajectories" . 

In order to compute the energy profile we write Hamil- 
tonian Q in terms of local energy density operators H n : 

H n = -Q<J n <7 n+1 + - ■ (ct„ + (T„+i) . (6) 

The local Hamiltonian H n (defined for < n < L — 2), 
gives the energy density between the n-th and (n + l)-th 
spins. In terms of eq. © the Hamiltonian of the system 
can be rewritten as 

W = £#n + ^i + ■ (7) 

n=0 

First we have performed equilibrium simulations in or- 
der to show that time averaged expectation values of the 
local energy density can be used as a consistent canon- 
ical local temperature. To this end we set the left and 



right reservoirs to the same temperature T. In Fig. [21 we 
plot the mean local energy E = (H n ) where the average 
is taken over time and over the L — 5 central values of 
the profile, for the chaotic and for the integrable chain 
(inset). For low T, E saturates to a constant which, to- 
gether with the energy profile (H n ), is determined by 
the ground state. However, for larger T > 1, the en- 
ergy profile is constant within numerical accuracy, and 
numerical simulations give E <~ — 1/T, all results being 
almost independent of L for L > 6. The numerical data 
for E(T) can be well approximated with a simple cal- 
culation of energy density for a two-spin chain (L = 2) 
in a canonical state at temperature T (smooth curves), 
namely E can (T) = tr H e- H °/ T /tr e- H °/ T . Therefore, if 
the temperatures of both reservoirs are in high T regime, 
then we can define the local temperature via the relation 
T oc —1/E. We stress that equilibrium numerical data 
shown are insensitive to the nature of dynamics (con- 
sistent with results of Ref . [l4^ ) . whether being chaotic, 
regular or intermediate. 

In Fig. we show the energy profile (H n ) for an out of 
equilibrium simulation of the chaotic chain. In all non- 
equilibrium simulations, the temperatures of the baths 
were set to T\ = 5 and T T = 50. After an appropriate 
scaling the profiles for different sizes collapse to the same 
curve. More interesting, in the bulk of the chain the 
energy profile is in very good approximation linear. In 
contrast, we show that in the case of the integrable (inset 
I) and intermediate (inset II) chains, no energy gradient 
is created. 
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FIG. 4: Size dependence of the energy current in the chaotic 
chain with 21 = 5 and T T = 50. We show J/ AE (open circles) 
and J /AT (open squares). The dashed lines corresponds to 
1/AL scaling. In the inset, the size dependence of the en- 
ergy current is shown for the integrable (solid circles) and the 
intermediate (solid squares). 
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FIG. 5: Infinite temperature current-current time correla- 
tion function £(i), for two chaotic chains of size L — 18 (tick 
dashed line) and L = 24 (thick solid line) and for two inte- 
grable chains of size L = 18 (dashed line) and L = 24 (solid 
line) . 



We now define the local current operators through the 
equation of continuity: dtH n — i[H, H n ] = — (J n +\ — J n ), 
requiring that J n = [H n ,H n -i]. Using eqs. JBJ and Q 
the local current operators are explicitly given by 




Jn = h x Q (a z r: 



-i) 



1 < n < L - 2. 



(8) 



In Fig. 0] we plot J/ AE as a function of the size L of 
the system for sizes up to L = 20. The mean current J 
is calculated as an average of (J n ) over time and L — 8 
central values of n. The energy difference was obtained 
from the energy profile as AE = (Hi,— 5) — (-H3). Three 



spins near each bath have been discarded in order to be 
in the bulk regime. Since AL = L — 8 is an effective size 
of the truncated system, the observed 1/AL dependence 
confirms that the transport is normal. Moreover, also 
the quantity J/ AT, where AT = -l/(H L - 5 ) + l/(H a ), 
shows the correct scaling with the size L. On the other 
hand, in integrable and intermediate chains we have ob- 
served that the average heat current does not depend on 
the size J oc L , clearly indicating a ballistic transport. 

As a consistency check, we have also studied the energy 
current-current time correlation £(f) = (J(t)J) of the 
isolated system for (3 — > which, for high temperature, 
gives the conductivity asK = (3L J* °° £(t)dt. In Fig. 0we 
show that for the chaotic chain exhibits fast, possi- 
bly exponential decay while for the intermediate case we 
observe slow decay, possibly to a finite plateau, giving 
rise to ballistic transport. In the integrable case is 
constant since J(t) is a conserved quantity. 

In this paper we have shown that Fourier law of heat 
conduction can be derived from the pure quantum dy- 
namical evolution without any additional statistical as- 
sumption. This possibility is strictly related to the onset 
of quantum chaotic behavior. 
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